sun <- read_rds("~/GoogleDrive/PhD/Paper Three/mob_multicity.rds") %>% filter(codigo_destino != codigo_residencia & dia == "Sunday") %>% na.omit()
wed <- read_rds("~/GoogleDrive/PhD/Paper Three/mob_multicity.rds") %>% filter(codigo_destino != codigo_residencia & dia == "Wednesday") %>% na.omit()
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = sun)
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = wed)
s1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = sun)
w1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = wed)
lc2_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))
s2_lc <- lmer(formula = lc2_formula, data = sun)
w2_lc <- lmer(formula = lc2_formula, data = wed)
s3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = sun)
w3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = wed)
all_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))
s4_all <- lmer(formula = all_formula, data = sun)
w4_all <- lmer(formula = all_formula, data = wed)
AIC(s_basic,s1_lc,s2_lc,s3_wea,s4_all)
## df AIC
## s_basic 7 552433.7
## s1_lc 9 552144.1
## s2_lc 19 537242.1
## s3_wea 10 551799.6
## s4_all 22 536568.5
AIC(w_basic,w1_lc,w2_lc,w3_wea,w4_all)
## df AIC
## w_basic 7 1144243
## w1_lc 9 1143452
## w2_lc 19 1129340
## w3_wea 10 1142632
## w4_all 22 1127414
summary(s4_all)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) +
## I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip +
## scale(uf_cont_destino) + scale(parks_destino) + scale(agriculture_destino) +
## scale(natural_destino) + scale(industry_destino) + scale(water_destino) +
## scale(uf_cont_residencia) + scale(parks_residencia) + scale(agriculture_residencia) +
## scale(natural_residencia) + scale(industry_residencia) +
## scale(water_residencia)
## Data: sun
##
## REML criterion at convergence: 536524.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2995 -0.6948 -0.0535 0.6546 4.3021
##
## Random effects:
## Groups Name Variance Std.Dev.
## city (Intercept) 0.1374 0.3707
## Residual 0.4590 0.6775
## Number of obs: 260426, groups: city, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.616e+01 1.985e-01 131.821
## scale(dens_destino) -1.873e-01 2.400e-03 -78.037
## scale(dens_residencia) 1.227e-02 2.379e-03 5.158
## log(dist) -4.536e+00 3.259e-02 -139.216
## I(log(dist)^2) 2.215e-01 2.068e-03 107.117
## temp 1.008e-03 8.335e-04 1.210
## I(temp^2) -1.162e-04 1.806e-05 -6.433
## precip -2.297e-04 2.346e-05 -9.790
## scale(uf_cont_destino) 2.054e-02 1.766e-03 11.629
## scale(parks_destino) 5.549e-03 1.655e-03 3.352
## scale(agriculture_destino) 9.500e-02 2.316e-03 41.024
## scale(natural_destino) 8.498e-02 2.467e-03 34.446
## scale(industry_destino) 1.683e-02 1.919e-03 8.766
## scale(water_destino) 9.044e-02 1.556e-03 58.132
## scale(uf_cont_residencia) -1.456e-02 1.773e-03 -8.213
## scale(parks_residencia) 3.121e-02 1.646e-03 18.954
## scale(agriculture_residencia) 6.953e-02 2.150e-03 32.346
## scale(natural_residencia) 8.835e-02 2.233e-03 39.564
## scale(industry_residencia) 1.325e-02 1.854e-03 7.146
## scale(water_residencia) 2.749e-02 1.508e-03 18.229
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
summary(w4_all)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) +
## I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip +
## scale(uf_cont_destino) + scale(parks_destino) + scale(agriculture_destino) +
## scale(natural_destino) + scale(industry_destino) + scale(water_destino) +
## scale(uf_cont_residencia) + scale(parks_residencia) + scale(agriculture_residencia) +
## scale(natural_residencia) + scale(industry_residencia) +
## scale(water_residencia)
## Data: wed
##
## REML criterion at convergence: 1127370
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6722 -0.7201 -0.0640 0.6798 4.3742
##
## Random effects:
## Groups Name Variance Std.Dev.
## city (Intercept) 0.2811 0.5302
## Residual 0.4435 0.6659
## Number of obs: 556646, groups: city, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.622e+01 2.353e-01 111.429
## scale(dens_destino) -2.262e-01 1.572e-03 -143.901
## scale(dens_residencia) 2.899e-02 1.596e-03 18.160
## log(dist) -4.681e+00 2.283e-02 -205.074
## I(log(dist)^2) 2.415e-01 1.410e-03 171.240
## temp 2.065e-02 6.044e-04 34.160
## I(temp^2) -5.218e-04 1.335e-05 -39.082
## precip 9.593e-05 1.791e-05 5.355
## scale(uf_cont_destino) -8.779e-04 1.218e-03 -0.720
## scale(parks_destino) -7.387e-02 1.127e-03 -65.567
## scale(agriculture_destino) 3.981e-03 1.502e-03 2.650
## scale(natural_destino) 1.343e-02 1.577e-03 8.518
## scale(industry_destino) 3.213e-02 1.251e-03 25.672
## scale(water_destino) 3.356e-02 1.007e-03 33.321
## scale(uf_cont_residencia) -1.010e-02 1.169e-03 -8.641
## scale(parks_residencia) 1.091e-02 1.129e-03 9.667
## scale(agriculture_residencia) 3.713e-02 1.447e-03 25.663
## scale(natural_residencia) 7.813e-02 1.514e-03 51.590
## scale(industry_residencia) 1.061e-02 1.231e-03 8.622
## scale(water_residencia) 1.562e-02 9.661e-04 16.163
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
s_complex <- s4_all
w_complex <- w4_all
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "B"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "B"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "B"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "B"))
sun_pred <- sun %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))
map <- st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp") %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
ggarrange(p1,p2,p3, nrow = 1)
Log of mean daily flows in both directions
ggarrange(p4,p5,p6, nrow = 1)
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "M"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "M"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "M"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "M"))
sun_pred <- sun %>% filter(city == "M") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "M") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))
map <- st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp") %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
ggarrange(p1,p2,p3, nrow = 1)
Log of mean daily flows in both directions
ggarrange(p4,p5,p6, nrow = 1)
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "V"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "V"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "V"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "V"))
sun_pred <- sun %>% filter(city == "V") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "V") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))
map <- st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp") %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
ggarrange(p1,p2,p3, nrow = 1)
Log of mean daily flows in both directions
ggarrange(p4,p5,p6, nrow = 1)
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "S"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "S"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "S"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "S"))
sun_pred <- sun %>% filter(city == "S") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "S") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))
map <- st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp") %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Land Cover Data Model")
ggarrange(p1,p2,p3, nrow = 1)
Log of mean daily flows in both directions
ggarrange(p4,p5,p6, nrow = 1)